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Abstract 

By combining the first unquenched lattice computation of the B-meson binding en- 
ergy and the two-loop contribution to the lattice HQET residual mass, we determine 
the MS 6-quark mass, mb{rnb)- The inclusion of the two-loop corrections is essential 
to extract mb{mb) with a precision of ©(Ag^-.^/m;,), which is the uncertainty due to 
the renormalon singularities in the perturbative series of the residual mass. Our best 
estimate is rnb(jnb) = (4.26 it 0.09) GeV, where we have combined the different errors 
in quadrature. A detailed discussion of the systematic errors contributing to the final 
number is presented. Our results have been obtained on a sample of 60 lattices of size 
24^ X 40 at /? = 5.6, using the Wilson action for light quarks and the lattice HQET for 
the b quark, at two values of the sea quark masses. The quark propagators have been 
computed using the unquenched links generated by the T^L Collaboration. 
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1 Introduction 



Quark masses are fundamental parameters of QCD that cannot be determined by theoretical 
considerations only and, due to the confinement of quarks inside hadrons, cannot be directly 
measured. Quark masses can, however, be introduced as short-distance effective couplings. 
As such, they are scale and scheme dependent quantities, the values of which depend on the 
adopted definition. Nonetheless, quark masses are very important for phenomenology since 
they enter many theoretical predictions of physical quantities such as CKM matrix elements, 
6-hadron inclusive semileptonic decays, total widths, etc. This is the reason why, in the last 
years, much effort has been devoted to accurately determine their values. 

It is useful to classify quarks into two classes: light quarks, the masses of which are lower 
or of the order of Aqcd (the u, d and s quarks are light), and heavy quarks, with masses 
larger than Aqcd ( b and t quarks are heavy and, to some extent, the charm quark c too). 
Light-quark masses are extracted from hadron spectroscopy using lattice QCD simulations [0 
and also QCD sum rules and r decay data Heavy-quark masses can be extracted from 
the properties of hadrons containing heavy quarks: the S-meson spectrum from the lattice 
HQET |, |, the T (or J/^p) spectrum with lattice NRQCD [§, § or QCD Sum Rules [§-^ 
and mass effects in 3-jets bbg events ||15|, |16|. 



In this paper, we present the first unquenched HQET lattice calculation of the b quark mass. 
The idea (see also [|17|) is to combine the HQET unquenched lattice computation of the 



S-meson binding energy [|I^ with the recent next-to-next-to-leading ( NNLO ) perturbative 
calculation of the matching of the MS quark mass to its lattice HQET counterpart [^]. We 
stress that both unquenched simulations and NNLO matching are necessary ingredients to 
improve the accuracy of previous results [Q]. The former is necessary to control potentially 
large vacuum-polarization contributions to the S-meson propagator. The latter is crucial to 
reduce renormalon ambiguities in the continuum- lattice matching . After a careful analisys 
of the systematics errors, our best result is 



nibimb) = (4.26 ± 0.03 ±0.05 ± 0.07) GeV 



where the first error is statistical; the second is the systematic error from the spread of 
values due to the use of different time intervals, fitting methods, smearing types and cube 
sizes for the interpolating operators, the dependence of the results on the mass of the sea 
quarks, the calibration of the lattice spacing and an evaluation of the l/rrib corrections; the 
third is an estimate of the error due to the uncertainties in the values of and to the effects 
of higher-order terms in eq. (piSl) . A detailed discussion of the different errors can be found 
below. 

The paper is organized as follows: in Sect. 2, we briefly describe our method and give the 
main formulae we used; in Sect. 3, we discuss the lattice computation of the binding energy; 
details of the simulation and numerical results are presented in Sect. 4, where we also discuss 
the procedure used for analyzing the unquenched lattice data; in Sect. 5, we carefully study 
the different sources of systematic errors in our results. Finally, in Sect. 6, we present our final 
numbers and compare them with other recent determinations. 
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2 The method. 



The key idea to determine the 6-quark mass consists in matching the propagator in QCD to 
its lattice HQET counterpart [0]. As shown below, this matching allows us to relate the pole 
mass to the binding energy and to the physical mass of the B meson. The renormalized MS 
6-quark mass at a given scale /i can then be obtained from the pole mass by using perturbation 
theory. In this section we briefly recall the formulae relevant to our study. 

Lattice HQET is an effective theory of QCD. The relation between the inverse 6-quark 
propagator in QCD, S~^, and its lattice HQET counterpart, Sj^^, can be written, to lowest 
order in l/nif,, as 

(^] Sp\p,m,-fi) ^ (^] (^] = C{fia,as) (^] Se\{v ■ k)a) , (1) 



where Sp^ is the projected 6-quark propagator, /i the renormalization point, a the lattice 
spacing, p = ruhV + k the momentum of the 6-quark, v its velocity and k the residual 
momentum. C(/ia, a^) is the relevant Wilson coefficient. It contains all the mass dependence 
of the right hand side of eq. ([1|) since, by construction, the HQET propagator is independent 
of the b quark mass. It should be noticed that in order for the HQET to be applicable, k must 
satisfy the condition -C m^. In writing eq. (0) we have chosen as expansion parameter the 
quark mass appearing in the original propagator, namely m^. 

The procedure to find C is well known: calculate the 6-quark propagator in QCD and in 
the lattice HQET to a given order in a^, expand the former in inverse powers of m;, to a given 
order (lowest order in our case), and finally compare both expressions at a fixed scale /i (with 
II ^ ^qcd) to extract C(/ia, a^). Renormalization group can then be used to evolve this 
function to any scale. 

To illustrate and clarify the key points of this strategy, we briefly sketch the derivation of 
our master formula for the 6-quark mass to 0(a<j) and then we extend our equation to include 
higher orders in perturbation theory. The inverse quark propagator in QCD can be written 
in the form 

iS~'^{p,mb;fi) = 'J) - rrib + Si(p^mb) + {'J) - m6)S2(p^mfe). (2) 

It is very easy to calculate the self-energy form factors Si and S2 to one loop in some renor- 
malization scheme and for a fixed gauge. By writing the 6-quark momentum in the B meson, 
p, as J9 = mf,f -|- k and expanding in powers of l/m^, one finds 



Sp^{p,mb; fi) = rrib 



m 



pole 



(vk) 



(3) 



where C2 is a scheme dependent constant (the expression of which is irrelevant for our dis- 
cussion) and 7m and lo-k are the scheme-independent one-loop anomalous dimensions of the 
mass and operator D ■ k, respectively. 

The pole mass, m^"'*^, is defined as the position of the pole of the propagator S^^, at a 
given order in perturbation theory, 

Sp^i:^ = ml°^^,mh]fi) = 0. (4) 
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To one loop, the explicit calculation of eq. (j^) gives 



pole 



rrib 



1 + i-fm In 



Cl 



(5) 



In order to implement the matching S*^^, the propagator of the lattice HQET, must be 
evaluated at the same order in perturbation theory 



Sl\{vk)a) = (vk) 



1 + asia) 'jD-k In 



-2{v ■ k)a^ 



do 



asia) 



(6) 



where d2 is a scheme dependent constant. From eq. (|]) we learn that an additive, linearly 
divergent mass term is generated on the lattice: the so-called residual mass, dm. Inserting 
eq. @ into eq. (^) and taking into account the expression of the pole mass in eq. (1), we 
obtain 

Sp^{p,mb] fi) = TJib - ml"^" + a^ifi) h C{^a,as) S£\{v ■ k)a) , (7) 



where the Wilson coefficient has the form 

C{fia,as) = 1 + as(yu) {in-k In (//a) + C2 - ^2) 



(8) 



and we have used the fact that the difference between ^^(/i) and as{a) is 0{af). Comparing 
eq. (|l|) with eq. (0), the important relation between the HQET expansion mass parameter 
and the pole 6-quark mass can be derived 



m' 



pole 



rrib + OisM — - = mb + 5m. 
a 



(9) 



To lowest order in l/rrib, the HQET mass formula can now be used to eliminate the un- 
known expansion parameter, rrib, by expressing it in terms of the physical mass of a fe-hadron, 
specifically the 5-meson, and the non-perturbative binding energy, S, which is independent 
of rrib, 

Mb = rrib + S + 0{l/mb) (10) 
Using the equation above, we get 



rri' 



pole 



Mb - £ + Sm + 0{l/mb) . 



Finally, the pole mass is converted into the MS mass through the well-known one-loop per- 
turbative relation 



mbirrib 



m 



pole 



1 - 



4 / asirribY 



TV 



X 

Mb - S + asijrib) — 
a 







l^sirrib)^ 




3 





Oil/rUb) . 



(12) 



We stress that mb{mb) is obtained from the non-perturbative quantity £ combined with the 
perturbative calculation of lattice (q;<jXo) and continuum (4/3 Q!s(rnfe)/7r) coefficients. 
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The generalization of eq. ([T^ ) to higher orders is straightforward. One gets 



pole 



1 + E 



n+1 



n=0 



TT 



+ ^(l/mfc) 



Ms - £ + ^(«,(mfe)) 

n=0 



— . ^^n+1 ^'-ra 

a 




n+1 



(13) 



which is the master equation of our analysis. -D„ and X„ are constants which depend on the 
number of flavours, rij, the masses of the active quarks and the lattice action used for the 
light quarks (see below). 

The procedure used to calculate fnbijnb) is the following [Q: 

• compute the binding energy of the B meson in lattice units, aS, using the HQET for 
the heavy quark and a given discretized action (Wilson, Alpha, Staggered) to describe 
the light-quark dynamics (we discuss in detail this computation in Sects. 3 and 4). 

• evaluate the value of the lattice spacing, a, from the light-hadron spectroscopy. 

• take the experimental value of the i?-meson mass, Mb, as input. 



• insert the values of these quantities in eq. (|13[) and obtain mf,{mf,) to a given order in 
perturbation theory and up to 0{l/mb) corrections. 

A few important remarks are in order at this point: 

1. The bare binding energy S is not a physical quantity since it diverges linearly as a — > 
and needs to be subtracted [0, |2^. Since both the pole and the MS mass are finite 
quantities, instead, the divergence of S is cancelled by the corresponding divergence of 
the residual mass 6m, expressed in terms of the constants Xn- 

2. In practice the cancellation is incomplete because we only know the values of few cons- 
tants Xn- Therefore, at a given order of the perturbative expansion, we cannot take the 
lattice spacing too small. 

3. The large- rij approximation shows that the perturbative series for 6m, and hence for 
the pole mass, suffers from renormalon singularities ||T7| , pOl . In other words, the coef- 
ficients Xn are expected to grow as const, x n! as ra oo. These singularities give 
rise to ambiguities of 0{Aqcd) in the sum of the perturbative series. A solution to this 
problem, which is the one adopted here, is to consider a short distance definition of the 
6-quark mass, such as the MS, frih, because it is free of renormalon ambiguities (up to 
0{AQfj^/mh) at the order at which we are working). 



In the expression of the MS mass, a delicate cancellation of renormalon singularities 
occurs: the renormalon of the series for 6m (with coeffcients Xn) is cancelled by the per- 
turbative expansion of the coefficient relating the pole and the MS mass (with coeffcients 

Dn). 
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5. In order to achieve the cancellation of renormalon singularities in eq. (|T3[), the same 
coupling constant has to be used in the expansion of 5m and in the relation between the 
pole mass and the MS mass [1T7| , For this reason, although we work at a fixed order 
of perturbation theory, we believe that the most reasonable choice is to expand both the 
continuum and lattice series using the same coupling constant. 

From the discussion above, it is clear that the precision of our results for the 6-quark mass 
at given a) is limited by the number of terms calculated in lattice {Xn) and in continuum 
{Dn) perturbation theory. The relation between the pole and the MS mass has been obtained 



to 0{al) by Gray et al. [0] and, recently, to 0{al) by refs. ||, |3 



Do 

D2 



4 

3' 



-11.6656 + 1.0414 J2 



^ rrij 



-157.116 + 23.8779n/ - 0.6527 ri 



(14) 



Note that the three-loop correction has been evaluated with massless quarks. 

As for the residual mass, 6m, it can be expressed in terms of the bare lattice coupling, oq, 

as 

6m 



oo Y' 
n=0 ^ 



The constant Xq is simply Xq - 

Xo = Cf 



-- Xo given by the three-dimensional integral 

1 r 1 



d-'k 



2 Eti sin^(A;./2) 



2.1173 



(15) 



(16) 



where Cp = {N"^ — 1)/2N and is the number of colors. Martinelli and Sachrajda have 
performed the calculation of Xi by extracting 6m from the exponential decrease of the expec- 
tation value of large Wilson loops with the perimeter [Q. More recently, Burgio et al. have 
obtained a preliminary estimate of X2 from large Wilson loops computed with the Numerical 
Stochastic Perturbation Theory (NSPT) on a 24^ lattice in the quenched approximation 
In summary, the results are 



^0 

X, 

X2 



2.1173 
11.152 + Uf 
73.5(9.2) 



0.282 + 0.035 C5H/ -0.391 4,^) 



;i7) 



where the value of the coefficient csw depends on the lattice fermions used in the simulation: 
for Wilson fermions csw = 0, for Clover-SW tree-level improved fermions csw = 1 and for the 
non-perturbatively improved ones, csw depens on j3 (see @] for details). The numerical value 
of X2 has been obtained in the quenched approximation (rij = 0) and thus it does not include 
fermion-loop effects. 

The next step is to express 6m in terms of the MS coupling as, i.e. to calculate the 
coefficients Xn from the Xn of eq. i\l7\) . The relation between as and ao can be written, to 
0{al), as: 



(/i) = tto + di{fia) al + d2{na) a^ + 



(18) 
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The pure gauge contributions to di and d2 have been calculated in ref. |25| and |26| respectively. 
The quark contribution to di for Wilson fermions can be found in and for improved 
fermions with generic csw in [H (see also ||2^). Unfortunately, the quark contribution to the 
two-loop coefficient d2 is still unknown. This calculation is necessary if X2 is used to obtain 
the N'^LO mass in the unquenched case. So far we have 

di{x) = ln(x) - 7;^ + 2.13573 iV 

Ztt l iV 

+ nj (-0.08413 + 0.0634 csvK - 0.3750 4iy) (19) 

17jV2 37r2 
d2{x) = di{xf - —— ln(x) + — — - 2.8626216 + 1.249116 iV^ 

IZTT^ 8iV^ 

where /3o = | (11 — 2n/-) and ^2 is given for rij = 0. By inverting eq. (|T^) for and 
inserting it in eq. (pTS]), for = 3 we get the values of the constants -^0,1,2^ 

Xo = 2.1173 

Xi = -1.30533 + 3.70677 In(mfea) 

+ nf{-0. 103872 - 0.224653 ln(m6a) - 0.0992368 csw + 0.402988 4^^) 

X2 = X2 + 6.48945 (-3.57877 + In(mba)) (3.29596 + In(mfea)) (20) 

Since the value of X2 is preliminary, we cannot really use it to obtain our final result. We 
will only show that, even in the quenched case, the result of ref. |^ is not precise enough to 
obtain a useful information. 



3 Lattice computation of S. 

The bare binding energy of the B meson, S, is measured on the lattice by studying the large- 
time behaviour of the B meson propagator The 6-quark is described by the discretized 
version of the HQET, 

Chqet = b{x) Dib{x) , (21) 
with the covariant derivative defined as 

D4 b{x) = U^{x) b{x + A) - b{x) , (22) 

where /t indicates the /x-direction and U^{x) is the link variable between the lattice sites x and 
x + [i. Light quarks, g, are simulated with some fermion action, in our case the Wilson action. 

It is well known that correlation functions involving heavy quarks suffer from a large 
contamination from higher-mass excitations to the lightest-state contribution, to which we 
are interested in. In order to improve the isolation of the lightest state, we use cube and 
double-cube smeared axial-current operators of size L^, as interpolating operators of the B 
meson [p9[| : 

A^(x) = 6(x) 75 g(x) , 

= EK^^)7/.75g(x), (23) 

i 



K^^) = II ^(^i)7M75g(2;j) 
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From the above operators we construct the two-point correlation functions 



CEfit) = X:(0|Af(f,t)Af t(0,0)|0), (24) 

X 

where R,R' stands for L (Local), S (Single smeared) and D (Double smeared) interpolating 
operators. The correlation functions are computed after rotating the links in the Coulomb 
gauge. 

At large time distances, C^^ (t) behaves as: 



CEfit) Zlzle-'\ (25) 



where 



By fitting the large time behaviour of C^^' (t) to eq. (p5|), the bare binding energy can be 
extracted. Due to contamination from excited states, the interpolating operators couple to 
the ground state in different ways so that the actual value of S has some dependence on the 
cube size and the smearing type used in the analysis. To obtain our best estimate, we compare 
different methods and account the mixing as a systematic effect contributing to the final error. 

We first use the Best Cube Method (BCM): we base our results on the best cube, defined as 
the one which yields the largest and flattest effective- mass plateau In practice, we proceed 
as the following example illustrates. Consider the double smeared interpolating operators. For 
t > train and for each L^, we fit C^^{t) and C^^(t) to a single-state propagator (^). t^m is 
taken as the time at which we start observing a plateau for both the effective mass ^E^^ and 
the ratio R?^^^^^ (we call BC DL/DD method the case where the smearing is in the sink) or, 

alternatively, for the effective mass AEff and the ratio R^f^^^ (we call BC LD/DD method 
the case where the smearing is in the origin). Effective masses and ratios are defined by 

ryRR' {+\ yR yR' 

r^RR'/NN' f.-, ^ ^Ls y'^) ^Ls ^Ls 

^La C^^'it) Z^' ' 

Ls ^ ' Ls Ls 

For t > tmin, the ground state is assumed to have been isolated. The next step is to combine the 
exponential fit for C2^{t) {C^^{t)) and the average value of the ratio 'R!l^^^^{t) {K^^^'^'^ [t]) 
in the plateau region to obtain both £ and Z^, the matrix element for the local axial current. 
Similarly, the method is applied to single smeared operators. 

A drawback of the BCM is that, in practice, we have only few different cubes at disposal 
(essentially only two cubes are really useful, Lg = 7 — 9, as suggested by earlier lattice studies). 
Nevertheless, this method is able to give a reasonable isolation of the lightest state and an 
accurate determination of the binding energy (although the method is less efficient for the 
determination of the matrix elements of the local axial current). 

In order to improve the accuracy of our analysis, we have also used the Multifit Method 
which consists in performing a global fit of the data for all cube sizes and smearing types by 



7 



minimizing the total 1 29]. Consider, for example, the double-smeared operators. In this 
case, Xtotai is defined by 



'X.total 





cE!'{t)DATA-cE!'{t) Y 



{2i 



where cr^f^^^^ is the jacknife error of the data points. We also impose the consistency condition 
that the binding energy S should be the same for all smearing types and cubes sizes. Moreover, 
in order to reduce the effect from excited states, it is convenient to fit the data by including 
the contribution of at least one excited state to C^^ {t). 



4 Strategy of the unquenched analysis. 

Before giving details on the simulation parameters, the calibration of the lattice spacing and 
the extraction of the binding energy S, we discuss the general strategy followed in the analysis 
of the unquenched data. This is a crucial issue due to the confusion which, we think, exists 
in the literature. 

Our main point, which is justified below, is that the correct procedure consists in per- 
forming independent quenched-like calculations of all quantities for each fixed value of kgea, 
including the spectroscopy, the calibration of the lattice spacing a and the calculation of the 
relevant matrix elements [|18[. Only when all the quantities are expressed in physical units. 



the results can be extrapolated in the sea quark mass. Thus, one ends up with a different set 
of lattice parameters for each kgea, such as the critical kappa, kcr, the light quark masses ku 
and kg, the lattice spacing and so on. Extrapolations in the valence quark mass, ky, should 
also be performed at fixed k^ea, without ever mixing up the extrapolation in the valence and 
sea quark masses, which must remain distinct steps of the analysis. In all respects, the value 
of ksea is an external "field" which controls the link dynamics. 

The argument is the following. A change in the value of the sea quark mass(es) modifies 
the value of the effective coupling constant, because the latter receives contributions from 
virtual-quark loops. A change of the coupling constant may induce a rapid variation of the 
value of the lattice spacing which depends exponentially on a^. Therefore, strictly speaking, 
lattice results for different kgea correspond to different lattice dynamics and are not directly 
comparable. Only when the results have been converted to physical units, by using the lattice 
spacing extracted for each kg^a, comparisons and extrapolations are possible. We stress again 
that a combined (in valence and sea quark masses) chiral extrapolation of lattice quantities, 
as for example the quark masses, may produce incorrect results because in this way we are 
mixing results corresponding to different values of the lattice spacing and all the parameters 
of the extrapolation do depend on kgea through a. 

Having explained our strategy, we turn to the numerical results for the binding energy. 
We have performed an unquenched computation of S with two degenerate sea quarks at two 
values of their mass, fc^ea = 0.1575 and kgea = 0.1580. The heavy and hght quark propagators 
have been computed using the set of unquenched link configurations generated by the T^L 
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Table 1: Simulation parameters of our run to extract the binding energy S. The values of 
the lattice spacing, a, and the critical (kcr), light (ku) and strange (kg) Wilson parameters are 
given separately for each sea quark mass (corresponding to ksea)(see text). 



Simulation Parameters 


= 5.6 


V = 24^ X 40 


nf = 2 Nconf = 60 


kv 


0.1560, 0.1570, 0.1575, 0.1580 


ksea = 0.1575 


k,,a = 0.1580 


= 2.51(6) GeV 


= 2.54(6) GeV 


kcr = 0.15927(5) 


kcr = 0.15887(4) 


ku = 0.15920(5) 


k^ = 0.15880(4) 


ks = 0.15747(12) 


ks = 0.15715(8) 



Collaboration. Details of the simulation can be found in ref. |^. Light quarks are simulated 
using the Wilson action whereas heavy quarks are described with the discretized HQET. The 
parameters of our run are given in Table |l]. The calibration of the lattice spacing has been 



performed using the K* — K lattice-plane method of ref. [^J • 

With the BC method, we find that the flattest and largest plateaus for AE^^ (t) and 

H^f (t) correspond to the cube size Ls = 7 and the correlation functions LD (smearing 
in the origin) and DD. With the Multifit Method, we also obtain that a two-state global 
fit of the correlations LD and DD describes very well the data. The agreement for single 
smeared interpolating operators is, instead, much worse. Therefore we base our results on 
double smeared operators since this is the most efficient way of isolating the lightest state. 
Our best estimates of the values of the binding energy S are 

aSs^ = {0.588(11)(5), 0.606(15)(2)} 

aSs, = {0.620(8)(5), 0.632(12)(2)} (29) 

for the two values of ksea = 0.1575 and 0.1580, respectively. The first error is statistical and 
the second systematic. The latter has been obtained from the spread of our results due to 
different time intervals of the fit, cube sizes and smearing types. Since a full account of the 



different methods and evaluation of the uncertainties can be found in ref. we do not give 
further details here. 



5 Sources of systematic error. 

Using eqs. (p!3|), (0), (^) and (p9|), we can readily obtain the value of the b quark mass. As 
the value of X2 is still preliminary and incomplete, we derive our results with the two-loop 
formula corresponding to the NNLO matching. In order to evaluate the systematic errors on 
these results, we carefully studied the different sources of uncertainties coming from the use 
of eq. (^) at this order: the value of as, higher-order perturbative corrections, input meson 
mass, method of extracting S, value of ksea and l/m^ corrections. In the following, unless 
stated otherwise, the central values in the tables correspond to kgea = 0.1580, S has been 
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Table 2: Dependence of mfy{miy) on the choice of Og. For Uf = 2 the value of Aqqd is 
explicitely given; the label Mz indicates that the physical value of as{Mz) has been used as 
input to compute agijnb) (see text). 



Dependence on 


Aqcd 


^/ 




fribimb) 


Mz 


5 


0.221 


4.31(3)(3) 


250 


2 


0.172 


4.25(3)(3) 


300 


2 


0.182 


4.26(3)(3) 


350 


2 


0.191 


4.27(3)(3) 



obtained with the Multifit Method, the input mass is the Bs meson mass and a hnear chiral 
extrapolation in to kg has been performed. In the tables, the first error is statistical and 
the second the systematic one based on the spread of the results due to the uncertainty in the 
lattice spacing, the different forms of writing eq. ( [I3| ) (see sect. 5.2) and the use of different 
time-intervals and fitting methods. 

5.1 Dependence of rnh{rnb) on a^. 

In order to obtain nibi^b) we have to choose the value of as to be used in the perturbative 



calculations. At all orders, to cancel the renormalon singularities in eq. (|13|), the same coupling 
constant has to be used for the lattice and continuum series. For this reason, although our 
calculation is truncated at we prefer to take the same coupling constant for both the 

lattice and the continuum cases. 

One possibility is to consider the physical value of obtained by running the experimental 
coupling as{Mz) = 0.118 down to as(mb) with Uf = 5. 

A second possibility is to account that our simulation has been performed with uj = 2 
and compute at the NLO with a (still to be determined) ^^cd- ^^e quenched case, the 
value of K'^cD ~ 250 MeV has been measured in ref. Since the physical value of Aqcd is 
expected to be larger than the quenched one, in the second case we have used the NLO value 
of as{mb) obtained by varying Aq^^ in the range [250, 350] MeV. 

The results for different values of as are presented in Table 0. The dependence on A!qcd 
is very weak: the maximum spread of the values is less than 20 MeV. The difference between 
the central values obtained with the Uf = 5 and Uf = 2 coupling is of about 50 MeV. We 
have taken this as a very conservative estimate of the error due to the choice of the coupling 
constant. 

Note that the value of as{mb) for A!qc£, = 300 MeV, as{mb) = 0.182, corresponds with a 
very good approximation to the arithmetic (and geometric) average of the quenched (n/ = 0) 
and physical {uj = 5) couplings (^^(mf,) = 0.15 and as{mb) = 0.22 respectively). For this 
reason our central value for rfibirnb) is that computed with ^^(mf,) = 0.182. 
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Table 3: Dependence of rnb{rnb) on higher orders in PT: (I) expanded and not expanded 
form of the master equation; (II) excluding or including preliminary (incomplete) NNLO 
contributions. Foruf = 2 the value of Kqcd is explicitely given; the label Mz indicates that 
the physical value of as{Mz) has been used as input to compute asijnb). 



Dependence on higher orders (I) 


^QCD 


^/ 


dsijrih) 


mh{mb) not expanded 


mb{mh) expanded 


Mz 


5 


0.221 


4.28(3)(1) 


4.34(3)(1) 


300 


2 


0.182 


4.24(3)(1) 


4.28(3)(1) 


Dependence on higher orders (II) 


^QCD 


^/ 


dsirflb) 


mb(mb) at 0{a'j.) 


mbirrib) at 0{al) 


Mz 


5 


0.221 


4.34(3)(1) 


4.16-4.34 


300 


2 


0.182 


4.28(3)(1) 


4.10-4.27 



5.2 Dependence of mb{mb) on higher orders. 

Eq. (|T^), which is used to evaluate fnj,(jnb) to C(a^) consists in the product of two factors. We 
can, then, organize the formula including (not expanding) or excluding (expanding to (9(a^)) 
the 0{al) terms arising from the product. We take the difference between these (formally 
equivalent) procedures, as an estimate of unknown higher-order terms in perturbation theory 
(PT). In Table ^, the values of mbijUb) obtained from the expanded and not expanded forms 
of eq. (p!3| ) are presented. For Uf = 2 the dependence of our results on higher orders is 30- 
40 MeV, for nf = 5 the difference is 60 MeV. From these spreads, we conclude that a fair 
estimate of the effect of expanding or not expanding eq. (|TB|) is ~ 50 MeV. From this estimate 
we assume from higher-order terms an error of ±25 MeV. 

As best estimate of rnb{mb), for each choice of the value of the coupling constant, we take 
the average of the results obtained with the not-expanded and expanded form of the master 
formula. In this way we have computed the central values given in tables I, I, I, I and 0. 

In order to get an independent estimate of the systematic uncertainty due to higher- 
order terms in the perturbative espansion, we also tried to compute mbirfib) using for X2 
the preliminary result of eq. (^). In the numerical calculations, obtained with the expanded 
form of eq. (|13]), we allowed X2 to vary in the la interval [64.3, 82.7], obtaining the range of 
masses given in Table |^ (only the central values are given). The values of Wibirfib) at order 
0{a1) are also given for comparison. There are huge numerical cancellations occurring in the 
calculation of X2 from X2 in eq. (^). For this reason, the difference between the NNLO and 
the (approximate) N^LO results varies from about zero to 180 MeV, depending on the value of 
X2- This quantity, even in the quenched case, is still affected by such a large uncertainty 
that it is impossible to use it for any realistic estimate. We urgently call for a more precise 
determination of X2 in both the unquenched and quenched cases. 

5.3 Dependence of miijUb) on the input B meson mass. 

Consistent values of mbirrib) should be obtained using as input either the or the Bg meson 
masses. The corresponding values of the binding energy, £b^ and £3^^ respectively, are used in 



11 



Table 4: Dependence of mi,{mb) on the input B mesons masses (see text). 



Dependence on Mb 


^QCD 




asirrih) 


mh{mij) from 5^ 


mi,{mh) from Bg 


Mz 


5 


0.221 


4.28(4)(3) 


4.31(3)(3) 


300 


2 


0.182 


4.23(4)(3) 


4.26(3)(3) 



Table 5: Dependence of mi,{mij) on the method to extracts (see text). 



Dependence on the method of extracting £ 


^QCD 


^/ 




LD/DD with = 7 


Multifit 


Mz 


5 


0.221 


4.34(2)(4) 


4.31(3)(3) 


300 


2 


0.182 


4.29(2)(3) 


4.26(3)(3) 



the two cases. This checks the lattice value of the Mb^-Mb^ mass splitting and the smoothness 
of our chiral extrapolation. Actually, for the Bs meson, the physical value of the strange quark 
mass, corresponding to A;^, is within the range of valence quark masses (see Table |^) and only 
a mild interpolation, rather than an extrapolation, is needed. For the Bet meson, instead, we 
have extrapolated almost to the chiral limit. In order to compute the pole mass we have taken 
Mb^ = 5.279 and Mb^ = 5.375 GeV [Q. In Table ^, we compare our results for the two cases. 
The results are nicely compatible and a small difference (~ 30 MeV) is observed in the two 
cases. 



5.4 Dependence of mb{mb) on the method for extracting S. 

In sect. 3, we discussed the two methods used to determine the binding energy S: the Best 
Cube method and the Multifit method. If the lightest state has been well isolated, both 
methods should give compatible results for the b quark mass. In Table ^, we present the value 
of rnb{rnh) for the BC method (obtained with LD/DD, smearing in the origin, and = 7, 
which gives the flattest and largest plateau) and for the Multifit method. Also in this case 
the results differ by ~ 30 MeV. 



5.5 Dependence of mb{mb) on ksea- 

As discussed before, in order to compare the values obtained for different ksea and attempt 
an extrapolation, we have first to convert the lattice quantities to physical units. In Table 
|, the values of mbijnf,) for either values of kgea are given. The dependence of our results on 
ksea, in the sea-quark mass region of our simulation, is small. Indeed, taking into account that 
they correspond to independent simulations, we are not able to observe any dependence on 
ksea within errors. Therefore, the only possible strategy is not to attempt an extrapolation 
in ksea and take the value at the lightest ksea, i-e. ksea = 0.1580, as the best estimate of the 
physical value of the b quark mass. The difference between the results at the two values of 
ksea is accounted as a systematic effect in the final error. 
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Table 6: Dependence of mb{mb) on kg^afsee text). 



Dependence on fc^ea 


^QCD 




asimC) 


mh^nib) with kg^a = 0.1575 


rriiyimb) with kgea = 0.1580 


Mzl 


5 


0.221 


4.34(2)(4) 


4.31(3)(3) 


300 


2 


0.182 


4.29(2)(4) 


4.26(3)(3) 



5.6 l/rrih correcctions to mb{mb). 

Our whole analysis is performed to the lowest order of the expansion in l/m^. This means 
that 0{l/mb) contributions to the relation between the QCD and the lattice HQET quantities 
in eqs. (|1]) and (p!0| ) have been neglected. We now make an estimate of the error introduced 
by these higher-order corrections. 

The HQET pseudoscalar mass formula including l/rUb corrections is given by 

MB = mb + S-^-^ + Oil/ml) , (30) 
2 nib 2 nib 

where the parameters Ai and A2 are matrix elements between B states of the kinetic and 
chromomagnetic operators 

{B\bD^b\B) 



Ai 



2M, 



B 



A2 . (ms-9B)m^ (3,) 



B 



with S the spin operator of the i?-meson and B the chromomagnetic field {Bi = | eijk Gjk)- 

It is straightforward to estimate A2 because which is related to the vector-pseudoscalar 
mass splitting, 



A2 = ^ (m|* - A'/|) ^ 0.12 GeV^ (32) 

The extraction of Ai, is, instead, more difficult as demonstrated by the spread of values 
obtained with different approaches: the lepton-energy spectrum in inclusive semileptonic B 
decays using Zero Recoil sum rules, QCD Sum rules, experiment data analysis and the HQET 
Virial Theorem (see and references therein). It has also been estimated on the lattice 



using the discretized HQET p, |3^ . Although biased by the lattice results, we prefer a small 



value for this parameter, in the absence of an accurate determination we let it to vary in the 
interval -0.5-0.0 GeV^. With this range, we find that the contribution of the l/nib corrections 
to the pole mass is at most ~ 30 MeV. Due to the theoretical uncertainties on Ai, we do not 
attempt to correct the l/m^ terms but include their effect as a sytematic error on the final 
result. 

5.7 Continuum limit of rnb{rnb) . 

To date, in quenched lattice simulations the binding energy aS has been computed at three 
values of P = 6.0, 6.2 and 6.4 |^ Q. In Table |^ we give the values of the 6-quark mass from 

^ Finite volume effects may be present in the results at = 6.4 since the lattice volume was rather small. 
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Table 7: Results formfy{rnh) at different values of the lattice spacing a. The first error is due 
to the uncertainties in the lattice determination of the binding energy, and the second is an 
estimate of the higher order perturbative corrections. 



Dependence on a in the quenched case 


f3 


a-\GeV) 




mb{mb) 


6.0 


2.0(2) 


0.61(1) 


4.34(5)(10) 


6.2 


2.9(3) 


0.52(1) 


4.29(7)(10) 


6.4 


3.8(3) 


0.460(7) 


4.25(7)(10) 



these quenched simulations for different lattice spacings. The quenched results are computed 
using the values of the binding energy from the APE Collaboration |2^ and NNLO quenched 
master formula with the coupling constant as{mb) = 0.15. The quenched values are very close 
to our new result with rtf = 2. Although one may argue that there is a (rather mild) tendency 
towards lower values as a decreases, with the present uncertainties we cannot attempt any 
extrapolation in a or realistic estimate of the discretization errors. 

6 Final result for nihijfih) and comparison with other de- 
terminations. 

We consider as best estimate of rfibijnb) the value obtained with the pole mass extracted by 
using the mass of the Bs meson and the binding energy £b^ measured on the lattice through the 
Multifit method at ksea = 0.1580, by averaging the results of the expanded and not expanded 
form of eq. (|l^) and by taking the NLO coupling constant computed at NLO with nj = 2 
and J^co = 300 MeV. Using the estimate of the different errors discussed in the previous 
section we then obtain: 

mb{rnb) = (4.26 ± 0.03 ±0.05 ± 0.07) GeV. (33) 

The first error is statistical. The second is the systematic error obtained from the spread of 
values due to the use of different time intervals, fitting methods, smearing types and cube sizes 
for the interpolating operators, the dependence of the results on the kgea, the calibration of 
the lattice spacing and the estimate of the l/rrib corrections. Finally, the third is an estimate 
of the error due to the uncertainties in the values of and to the effects of higher-order terms 
in eq. (|13D. We find that the latter is the most important source of error in the final result. 
For this reason a big effort must be done to compute the unknown N^LO contributions to 
the residual mass on the lattice and the NNLO matching coefficient between the lattice and 
continuum ctg in the unquenched case. On the numerical side, a non-perturbative calculation 
of ^QCD is ^ls° important. 

Our new result (^) modifies and improve the previous one obtained from quenched lattice 
simulations with NLO matching only 

fribimb) = (4.15 ± 0.05 ±0.20) GeV (34) 
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where the first error is due to the lattice systematics and the second is an estimate of higher 
orders. 

It is interesting to compare eq. (|33|) with recent values obtained with completely different 
approaches as mass effects in 3-jets bbg events, bb production cross-section and T spectroscopy. 
Our final result is in good agreement with most of NNLO estimates, as shown in fig.|l[ In the 
figure we also give our world average and error. This average has been obtained by using only 
the most recent NNLO determinations from T spectroscopy and lattice QCD, i.e. we did not 
use the results of refs. [0, |, ||, |^, either because they have been superseeded by more accurate 
calculations or because they are only computed at the NLO accuracy. The average is 

mb(mfe) = 4.23 ± 0.07 GeV (35) 

which corresponds to a relative error of less than 2% comparable to the precision on the top 
quark mass. The masses of the quarks of the heaviest and last discovered generation are, and 
will probably remain, the most accurately determined quark masses. 
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Determinations Year > 1994 



This work 
Gimenez 97 [3] 
Martinelli 99 [4] 

Davies 94 [5] 
AN Khan 99 [6] 
Narison 95 [7] 

Jamin 99 [8] 
Voloshin 95 [9] 
Melnikov 99 [10] 
Hoang 99 [11] 
Beneke 99 [12] 
Pineda 98 [13] 

Penin 99 [14] 
Rodrigo 99 [15] 
DELPHI 98 [16] 

AVERAGE 



I • 1 




4.0 



4.25 



4.50 



in,(mj,) GeV 



Figure 1: Determinations of mb{mb) from different methods compared with our result. Only 
references after 1994 have been included. 
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